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1.0  Introduction 

For  more  than  two  decades,  computers  have  been  used  to  differentiate 
functions  by  applying  the  basic  rules  of  calculus  to  the  symbols  appearing  In 
the  functions  [1],  Symbolic  or  analytic  differentiation  is  of  interest  because 
*  it  appears  to  be  the  most  efficient  approach  to  evaluate  the  derivative  of  a 

function  at  a  series  of  points.  The  symbolic  expression  for  the  derivative  is 
found  once,  and  that  expression  is  then  evaluated  many  times  with  little  effort. 
Numeric  differencing,  the  contrasting  numerical  method,  requires  that  at  each 
point  the  function  be  evaluated  at  least  twice,  and  the  derivative  is  then 
approximated  by  taking  the  slope  of  the  line  joining  those  points. 

In  1962  Hanson,  Caviness,  and  Joseph  [2]  surveyed  the  methods  that  had  been 
published  for  symbolic  differentiation  and  proposed  a  routine  for  analytic 
differentiation  which  included  a  wide  class  of  allowable  functions.  By  applying 
the  transformation  algorithm  of  Erchov  [3],  their  routine  eliminated  the  highly 
structured  input  requirements  of  earlier  methods.  The  major  detractor  of  Hanson's 
routine  is  that  the  first  and  second  derivatives  that  the  method  produces  are 
inefficient  because  derivatives  tend  to  increase  in  complexity  at  an  exponential 
rate. 

The  method  of  symbolic  differentiation  reported  here  retains  the  data  struc¬ 
ture  of  the  ordered  triple  proposed  by  Hanson,  et.  al.  as  its  skeleton,  but 
introduces  improved  algorithms  for  the  production  and  manipulation  of  the  ele¬ 
ments  of  that  data  structure.  Implementation  of  this  method  has  resulted  in  a 
routine  which  produces  "optimized"  code  for  partial  derivatives  comparable  to 
special  purpose  subprograms  written  by  users  for  evaluation  of  first  and  second 
partial  derivatives.  An  experimental  computer  program  has  been  developed  on 
which  computational  results  are  reported.  The  routine  has  also  been  interfaced 


with  a  nonlinear  programming  code  and  experimental  results  indicate  that  this 
method  is  a  desirable  alternative  to  other  methods  of  evaluating  first  and  second 
partial  derivatives  of  user-defined  functions. 

1.1  Motivation 

The  authors'  philosophy  is  that  computers  are  intended  to  aid  people  in 
the  solutions  of  their  problems,  not  create  new  ones.  An  implication  of  this 
philosophy  is  that  computer  program  execution  time  alone  is  inadequate  for 
evaluating  the  cost  of  problem  solutions.  Time  spent  by  humans  in  problem 
formulation,  data  preparation,  and  so  on,  frequently  provides  the  major  con¬ 
tribution  to  the  costs  of  problem  solution,  regardless  of  whether  such  costs 
are  itemized  in  dollar  terms.  Furthermore,  the  emphasis  on  and  attention  paid 
to  such  "human"  aspects  of  scientific  computation  should  and  must  increase 
heavily;  the  pace  of  solid-state  technology  is  accelerating  the  increase  in 
speed  and  capacity  of  computer  hardware  faster  than  the  hardware  costs  are 
increasing.  The  human  aspects,  whose  cost  trend  is  just  the  opposite,  will 
inevitably  shoulder  aside  hardware  considerations  as  the  prime  cost  and  time 
bottle-neck  in  obtaining  problem  solutions. 

The  development  of  the  symbolic  differentiation  routine  reported  here  was 
motivated  by  the  above  philosophy  and  the  fact  that  the  authors  were  experiencing 
substantial  time  delays  in  using  many  codes  for  solving  nonlinear  programming 
(NLP)  problems,  since  most  solution  methods  require  that  a  set  of  functions  and 
their  first  and  second  partial  derivatives  be  evaluated  at  arbitrary  points. 

Such  NLP  computer  programs  require  that  the  user  provide  subprograms  which  are 
able  to  evaluate  the  functions  (often  having  many  variables)  and  their  partial 
derivatives  at  any  point;  consequently,  a  user  is  compelled  to  know  a  programming 
language  such  as  FORTRAN,  or  know  someone  who  does,  in  order  to  prepare  a  problem 


for  solution.  Substantial  time  and  effort  are  also  required  to  code  the 
needed  routines  and  debug  them  once  written.  Even  if  standard  numerical  dif¬ 
ferencing  techniques  for  first  and  second  partials  are  used,  considerable 
Instability  may  occur  depending  on  the  interval  size  used  and  the  function 
being  evaluated  [4,  5].  The  purpose  of  this  paper  is  to  describe  procedures 
used  in  developing  standard  routines  which  will  eliminate  these  problems. 
Collectively  these  procedures  will  be  referred  to  as  a  symbolic  differentiation 
routine,  SYDIF. 

The  results  of  the  symbolic  approach  to  differentiation  and  its  computer 
implementation  are  not  restricted  to  nonlinear  programming  problems.  Any  com¬ 
puter  procedure  which  requires  that  users  provide  subprograms  for  evaluating 
function  and  their  partial  derivatives  can  be  simplified  by  the  symbolic  dif¬ 
ferentiation  approach.  Some  areas  of  application  include  graphical  interactive 
routines  which  could  be  used  to  aid  courses  in  calculus  by  displaying  iso-curves 
of  functions  and  their  partial  derivatives.  The  feasibility  of  using  SYDIF  in 
conjunction  with  introductory  calculus  courses  taught  in  a  Computer-Assisted 
Instruction  mode,  is  also  under  consideration,  since  students'  answers  could 
be  checked  exactly  by  matching  their  expressions  against  the  symbolic  derivative 
produced.  Other  possible  applications  include  areas  such  as  numerical  integra¬ 
tion  of  systems  of  differential  equations,  solution  of  nonlinear  equations  by 
Newton's  method,  unconstrained  optimization,  guidance  and  trajectory  systems, 
and  applications  of  linear  programming  that  iteratively  solve  non-linear  problems. 

2.0  Methodology 

The  symbolic  differentiation  algorithm  follows  the  outline  of  such  algorithms 
defined  by  Hanson  [2]  as  "In  step  one,  the  mathematical  expression  is  analyzed 
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and  the  order  in  which  the  appropriate  rules  of  elementary  differentiation  are 
to  be  applied  is  determined.  Step  two  consists  of  the  actual  application  of 
those  rules." 

2. 1  Data  Structures 

To  accomplish  step  one,  SYDIF  transforms  the  input  expression  into  an 
equivalent  parenthesis-free  expression  represented  by  a  tree-structure  speci¬ 
fying  the  order  of  the  operations.  This  data  stucture  was  chosen  because  it 
can  be  produced  efficiently,  evaluated  easily,  differentiated  in  a  straight 
forward  manner  using  the  basic  rules  of  calculus,  and  the  resulting  derivative 
can  be  optimized. 

Each  "branch"  of  the  data  structure  represents  a  subexpression  of  the  user's 
function  as  a  binary  operator  and  its  two  operands.  The  expression  (x^  +  X2)X^-  3 
is  represented  by  the  following  tree,  where  t  represents  exponentiation: 


where  (1)  in  line  2  refers  to  the  subexpression  in  the  previous  line  1.  By 
sequentially  evaluating  the  subexpressions  beginning  with  the  first  one,  the 
expression  can  be  evaluated.  The  value  at  the  root,  R  is  the  value  of  the 
expression.  The  above  data  structure  is  also  equivalent  to  Polish  (parenthesis- 
free)  notation  [7]  which  can  be  produced  efficiently  by  parsing  the  original 
infix  (normal  algebraic)  notation  expression  using  the  algorithm  described  below. 

2.2  Parsing  Algorithm 

The  parsing  algorithm  differs  from  the  usual  inf ix-to-postfix  transformation 
algorithm  in  its  manner  of  using  an  intermediate  push-down  stack.  The  parsing 
algorithm  causes  the  input  string  to  be  scanned  from  left  to  right.  Operands 
encountered  are  placed  on  an  intermediate  push-down  stack.  This  is  different 
from  algorithms  which  only  produce  Polish  notation  since  they  move  operands 
directly  to  the  output  string,  but  it  is  required  since  only  complete  subex¬ 
pressions  can  be  removed  from  the  stack.  The  algorithm  of  Hanson,  et.  al. 
produces  a  similar  data  structure,  but  achieves  it  by  using  Erchov's  right-to- 
left  scan;  consequently  the  entire  expression  must  be  input  before  processing. 
Furthermore,  more  elements  must  be  kept  in  the  stack  since  operators  of  equal 
precedence  must  be  put  onto  the  stack,  hence  a  larger  stack  is  required  and 
more  time  is  spent  clearing  the  stack  using  Hanson's  algorithm.  The  algorithm 
presented  in  this  paper  improves  on  the  previous  algorithm  by  using  a  left-to- 
right  scan  which  takes  advantage  of  the  normal  lef t-to-right  execution  of  opera¬ 
tors  of  equal  precedence  to  achieve  more  efficient  storage  management. 

The  SYDIF  parsing  algorithm  is  quite  simple  and  can  be  described  as  follows. 
Operands  are  placed  directly  in  the  intermediate  stack.  Operators  are  compared 
using  normal  FORTRAN  precedence  with  the  operator  nearest  the  top  of  the  inter¬ 
mediate  stack.  If  the  operator  from  the  input  string  has  precedence  greater 
than  the  precedence  of  the  operator  in  the  stack,  the  input  operator  is  put 
on  top  of  the  stack.  The  "("  always  goes  to  the  stack.  If  the  precedence  of 
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the  input  operator  is  Jess  or  equal  to  that  of  the  stack  operator,  the  top 
three  entries  are  moved  from  the  stack  to  the  table  of  triples,  with  the  top 
operand  moving  to  the  right  operand  position.  A  new  operand,  representing  the 
branch  of  the  tree  where  the  triple  was  placed,  is  put  on  the  top  of  the  stack. 
The  input  operator  which  caused  the  ejection  from  the  stack  i9  then  compared 
with  the  top  operator  remaining  on  the  stack,  and  the  process  repeats  until 
the  end  of  the  expression  is  encountered.  The  end-of-expression  symbol  is  of 
equal  precedence  to  the  begin-expression  symbol,  both  the  lowest  possible,  hence 
the  stack  is  cleared.  Errors  in  the  syntax  of  the  expression  are  easily  recog¬ 
nized  by  the  parsing  algorithm. 

2.3  Differentiation  Algorithm 

The  expression  can  be  differentiated  after  being  transformed  into  a  tree- 
structure  because  such  a  data  structure  specifies  the  order  in  which  the  opera¬ 
tions  are  to  be  performed.  The  expression  is  differentiated  in  the  same  order 
that  it  is  evaluated;  beginning  with  the  first  triple  and  proceeding  sequentially 
until  the  root  is  reached.  The  new  "branches"  of  the  tree  representing  the 
derivatives  of  the  expression  are  attached  after  the  root  of  the  tree. 

Subexpressions  to  be  differentiated  are  classified  into  categories  according 
to  the  types  of  the  operands  and  operator.  Each  category  of  subexpression  is 
treated  as  if  it  were  an  assembly  language  macro-call,  the  macro  expansion  being 
the  partial  derivative  for  that  category  of  subexpressions.  The  parameters  of 
the  macro-call  are  the  triple  being  differentiated,  its  operands  and  operator, 
the  variable  of  differentiation,  and  pointers  to  the  derivative  of  any  reference 


to  preceeding  subexpressions.  Such  pointers  are  required  for  derivatives  of 
some  types  of  subexpressions  because  the  chain  rule  of  calculus  requires  that 


the  derivative  of  any  line  containing  a  reference  to  a  preceeding  subexpression, 
use  the  derivative  of  that  subexpression  in  obtaining  its  own  derivative. 

This  algorithm  improves  on  the  algorithm  of  Hanson,  et.  al.  by  recording 
the  actual  derivative  instead  of  a  pointer  to  the  derivative,  if  that  derivative 
is  a  single  operand.  By  later  referencing  the  actual  derivative,  approximately 
ten  percent  of  the  subexpressions  produced  for  second  partial  derivatives  are 
eliminated.  Such  a  recording  procedure  also  is  a  prerequisite  for  further 
optimization. 

2.4  Optimization  Algorithm 

The  underlying  purpose  of  the  symbolic  differentiation  routine  is  the 
production  of  efficient  code  for  the  evaluation  of  partial  derivatives.  Optimi¬ 
zation  actually  starts  with  treating  each  line  as  a  macro  call,  resulting  in  as 
few  lines  being  created  as  possible.  Later  specific  optimization  techniques 
are  applied,  both  as  triples  are  being  added  to  the  tree,  and  again  after  the 
entire  expression  or  partial  derivative  is  created  as  described  below. 

As  each  line  is  added  the  transitive  and  commutative  properties  of  the 
pairs  of  operators  (+,  -)  and  (*,  /)  are  used  to  fold  subexpressions  such  as 
C  *(R/C^)  into  *R2«  (where  and  are  constants  and  R  is  a  reference  to 
a  variable  or  previous  triple  and  =  C^/ is  a  new  constant) .  *R  has  a 
simpler  derivative  than  C^*(R/C2)  because  the  latter  requires  differentiation 
of  R/C2«  Thus  folding  shortens  both  the  expression  being  created  and  later 
partial  derivatives. 

Redundant  lines  are  removed  after  the  entire  expression  or  partial  deriva¬ 
tive  is  created.  The  SYDIF  optimization  algorithm  compares  each  triple  against 
previous  triples;  if  a  match  is  found,  the  duplicate  is  removed  and  all  later 
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references  to  the  removed  triple  are  changed  to  point  to  the  first  occurence  of 
the  triple.  The  comparison  process  starts  with  the  second  triple  being  checked 
against  the  first  and  is  completed  when  the  subexpression  at  the  root  has  been 
compared  against  all  the  preceeding  ones.  Triples  which  have  operators  +  and  * 
are  commuted  for  the  test;  special  checks  are  made  for  cases  such  as  (R  +  R)  = 

2  *  R  or  (0  -  R)  =  -1  *  R  by  also  checking  for  the  alternate  form. 

Dead  code  elimination,  (the  final  phase  of  optimization)  consists  of 
purging  superfluous  subexpressions  accumulated  in  creating  expressions  and  their 
partials.  The  algorithm  recursively  eliminates  all  triples  not  referenced 
later,  starting  at  the  root  of  the  tree.  In  the  above  example,  once  C^*R  is 
created,  there  may  be  no  later  references  to  the  triple  R/C2,  implying  that 
the  operation  need  not  be  present.  Dead  code  elimination  would  remove  the 
triple  R/C2  from  the  tree. 

3.0  Experimental  Program  Implementation 

An  experimental  implementation  of  SYDIF  has  been  coded  in  FORTRAN  and  is 
being  tested  on  the  CDC  6600  at  The  University  of  Texas  at  Austin.  This  imple¬ 
mentation  was  written  to  interface  with  the  Sequential  Unconstrained  Minimiza¬ 
tion  Technique  (SUMT)  of  Fiacco  and  McCormick  [6],  but  is  equally  applicable 
to  any  problem  requiring  partial  derivatives  or  user-defined  functions.  The 
user  of  this  routine  is  ultimately  interested  in  two  aspects  of  the  program: 
how  easy  is  it  to  use  and  how  well  does  it  work? 

The  authors'  feel  that  users  want  a  familiar,  simple  input  form.  The  SYDIF 
implementation  allows  users  to  input  their  functions  in  FORTRAN  notation  with 
FORTRAN  operators.  Variables  are  XI  through  X100.  Real  constants  in  E-notation 
are  allowed,  as  are  integers.  Implied  multiplication  is  allowed  by  inserting 


the  multiplication  operator  (*)  for  missing  operators.  Many  other  options  were 


also  Included  in  order  Co  provide  users  with  a  natural  means  for  expressing 


functions. 

The  experimental  program  outputs  optimized  pseudo-machine  code  subprograms, 
equivalent  to  symbolic  partials,  for  all  first  and  second  partial  derivatives 
for  each  specified  function.  These  representations  can  then  be  evaluated  by 
interpretation  or  they  can  be  translated  directly  into  machine  code  for  direct 
execution  as  a  subprogram  by  a  particular  computer  (see  [7]  pp.  342-347).  The 
accuracy  of  the  partial  derivatives  produced  is  limited  only  by  the  accuracy 
of  the  computer's  floating  point  hardware  and  system  functions.  Computational 
experience  with  the  experimental  implementation  of  this  code  has  exceeded  the 
original  high  expectations.  The  optimization  not  only  produces  more  efficient 
code  for  evaluation;  but  surprisingly,  it  also  results  in  a  faster  differentia¬ 
tion  phase  by  eliminating  many  subexpressions  of  the  original  expression  and 
its  first  partial  derivative  that  otherwise  would  have  required  differentiation. 
The  optimization  also  reduces  data  structure  storage  requirements  (by  over  80% 
in  some  cases)  that  hamper  Hanson's  symbolic  approach  to  differentiation,  e.g., 
with  no  optimization,  Hanson's  algorithm  for  the  differentiation  of  the  expres¬ 
sion  (X^  +  X^)**X^  gives  132  subexpressions  for  all  first  and  second  partial 
derivatives  (see  Appendix  A  for  this  example)  .  SYDIF  with  its  optimization 
left  only  24  subexpressions  to  be  evaluated  (see  Appendix  B) . 

The  experimental  SYDIF  program  takes  40,000  words  of  memory  on  the  CDC  6600. 
Tust  runs  on  a  problem  having  thirty  functions  (of  varying  complexity)  took  only 
two  CPU  seconds  of  processing  to  generate  all  first  and  second  partial  deriva¬ 
tives.  Currently  a  machine- independent  FORTRAN  version  is  being  developed  and 
should  be  available  soon. 
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4 . 0  Application  of  SYDIF  in  Non-Linear  Mathematical  Programming. 

An  enormous  body  of  theoretical  and  computational  research  has  been  directed 
at  achieving  efficient  solution  procedures  for  nonlinear  programming  problems  of 
the  form: 

Minimize  F(x) 
subject  to: 

gi(x)  £  0,  i  =  1,  ....  I. 

and  h.  (x)  ~  0 ,  k  1 ,  . . . ,  K  • 
k 

Particularly  powerful  results  have  been  obtained  for  the  case  in  which  the 
problem  functions  are  once-or  twice-dif f erentiable  and  where  the  function  F(x) 
and  the  set  of  points  X  satisfying  the  constraints  are  convex.  This  paper  is 
not  intended  to  provide  a  bibliography  for  nonlinear  programming;  see  [8]  for 
references.  Chapter  8  of  [6]  discusses  the  Sequential  Unconstrained  Minimiza¬ 
tion  Technique  (SUMT)  algorithm  of  Fiacco  and  McCormick,  which  is  a  procedure 
both  elegant  and  known  to  be  effective  from  a  computational  viewpoint.  A  com¬ 
puter  program  (in  FORTRAN)  implementing  the  SUMT  algorithm  is  available  from 
the  SHARE  Library  under  distribution  number  SDA  3189  (see  also  [9]). 

As  can  be  seen  from  (1),  the  statement  of  a  general  nonlinear  programming 
problem  requires  one  to  specify  the  objective  and  constraint  functions;  there 
are  I  +  K  +  1  such  functions.  In  addition,  many  of  the  procedures  for  solving 
(1),  (e.g.,  SUMT,  Zoutendijk,  Newton,  Davidon-Fletcher-Powell,  etc.)  use  partial 
derivatives  of  these  functions  with  respect  to  each  of  the  variables  as  well  as 
the  Hessian  (matrix  of  all  possible  mixed  second  partials)  of  each  function; 
even  the  simplest  direct  search  methods  require  at  least  the  ability  to  calculate 
the  values  of  the  objective  and  constraint  functions  at  arbitrary  points. 

SYDIF  provides  a  user-oriented  input  capability  and  creates  the  tree  struc¬ 


ture  for  evaluating  the  original  expression  and  all  the  partial  derivatives. 
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The  applications  program  "calls"  evaluation  programs,  giving  the  routine  a 
vector  of  values  for  the  variables  and  which  function  or  which  partials  are  to 
be  calculated.  The  triples  for  that  expression  or  derivative  are  evaluated  one 
by  one,  starting  with  the  first  one.  As  the  value  of  each  triple  is  calculated, 
that  value  is  stored  in  an  auxilary  table.  The  entire  process  is  simple  and 
fast.  Running  in  such  an  interpretive  mode,  the  interfaced  version  of  SUMT 
took  only  twice  as  long  as  when  subroutines  for  partials  derivative  evaluation 
were  coded  by  hand.  The  turn  around  time  and  the  human  interaction  were  shortened 
as  expected;  SYDIF  took  less  than  one  tenth  of  the  user's  time. 

Actual  machine  code  for  a  particular  computer  could  have  been  created,  but 
the  authors  feel  that  the  portability  of  the  machine- independent  approach  is 
more  important  than  complete  efficiency,  especially  in  the  development  stage. 

Later  versions  producing  direct  machine  code  will  produce  code  that  takes  less 
execution  time  than  is  possible  even  with  user-written  routines. 

5.0  Summary 

It  was  necessary  to  look  at  the  area  of  creating  routines  for  first  and 
second  partial  derivatives  because  of  the  lack  of  a  user-oriented  interface 
between  mathematical  programming  codes  requiring  derivatives  and  its  users. 
Symbolic  differentiation  was  chosen  because  of  the  instability/computer  problems 
Inherent  in  numerical  differentiation  techniques  and  the  inability  of  numerical 
methods  to  overcome  the  human  problem  of  providing  the  user  with  a  desirable 
input  format  for  his  problem  statement. 

Having  chosen  the  symbolic  approach,  parsing,  differentiation,  and  code 
optimization  procedures  were  devised  in  order  to  improve  upon  the  computer  time 
and  storage  requirements  of  early  symbolic  approaches  [2].  These  new  procedures 
were  then  computational  compared  against  early  approaches  and  found  to  be  sub¬ 
stantially  better  in  both  time  and  space  requirements. 
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A  user-oriented  front-end  was  added  to  the  experimental  symbolic  differen¬ 
tiation  routine.  Next  this  routine  was  interfaced  with  SUMT  in  an  interpretive 
mode  and  compared  against  SUMT  solution  tiuie  with  hand-coded  evaluation  subroutines. 
Running  in  an  interpretive  mode  the  SUMT  version  took  only  twice  the  computer 
time  and  reduced  the  user  Interaction  time  90%.  If  the  user  finds  the  increase 
in  execution  time  to  be  prohibitive,  a  non- interpretive  (or  direct  executing) 
version  could  easily  be  written  for  a  particular  computer. 
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Number  of  lines 
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®lx5 
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5 

^F 

Ik 

X2X2 

a2? 

23 

ax2x5 

JL- 

8 

ax5 

a2F 

20 

*5*5 

TOTAL 

132  line 

required 
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Appendix  B 

derivative  code  with  optimization  (All  partials  shown) 


(1) 

X1 

+  X? 

(2) 

h 

**  X, 

*) 

(3) 

S 

-  1 

w 

L1 

**  L3 

CO 

LJ* 

X5 

(<0 

s 

+  (-2.0) 

(7) 

L1 

**  L,- 

6 

(8) 

L7 

*S 

(9) 

X5 

*  h 

(10) 

LOG  (L  ) 

(11) 

Ln*h 

(12) 

x5  *  L11 

(13) 

h  +  L12 

(1M 

X,.  -  1. 

(15) 

L1  **  LlU 

(16) 

L15  *  Xt> 

(17) 

L7  *  L1U 

(18) 

X5  *  L17 

(19) 

L10  *  L15 

(20) 

X5  *  L19 

(21) 

L15  +  L20 

(22) 

Log  (Lx) 

(23) 

L2  *  L22 

(2*) 

L22  *  L23 

JL. 

ax. 


x5(x1  +  x2)x5  '  1 


axlxl 


^XlX2 


(X1  +  x2)x5  '2(x5-i)x^ 


value  of  Lr 


*1*5 


—  -  Log  (Xx  +  X2)  (Xx  +  Xg)V 

X5  +  Log  (Xx+  Xg)X5-l 


■  X5(X1  +  X2)V 


2 

^F 

ax2x2 


&x2x5 


<X1  +  X2>  ‘  <X1  +  X2>A5 


Log  (Xx  +  Xg)  Log  (X]+  Xg)  (XX+  X2)X9 


&x5x5 


